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Motivated by applications in systems biology, we seek a probabilistic frame- 
work based on Markov processes to represent intracellular processes. We re- 
view the formal relationships between different stochastic models referred to 
in the systems biology literature. As part of this review, we present a novel 
derivation of the differential Chapman-Kolmogorov equation for a general mul- 
tidimensional Markov process made up of both continuous and jump processes. 
We start with the definition of a time-derivative for a probability density but 
place no restrictions on the probability distribution, in particular, we do not 
assume it to be confined to a region that has a surface (on which the prob- 
ability is zero). In our derivation, the master equation gives the jump part 
of the Markov process while the Fokker-Planck equation gives the continuous 
part. We thereby sketch a "family tree" for stochastic models in systems bi- 
ology, providing explicit derivations of their formal relationship and clarifying 
assumptions involved. 
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' 1 Introduction 

X" 

Systems biology is a merger of systems theory with molecular and cell biology. The key 
distinguishing feature of a systems biology approach is the description of cell functions 
(e.g. cell differentiation, proliferation, apoptosis) as dynamic processes [???]. There are 
two dominant paradigms used in mathematical modelling of biochemical reaction networks 
(pathways) in systems biology: the "deterministic approach", using numerical simulations of 
nonlinear ordinary differential equations (inch mass action type, power law or Michaelis- 
Menten models), and the stochastic approach based on master equation and stochastic 
simulations. 

Key references in the area of stochastic modelling are the books by ? ], ? ], ? ] and 
? ]. Most stochastic models presented in these references are derived on the basis of the 
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Chapman-Kolmogorov equation (CKE), a consistency condition on Markov processes, in 
the form of a system of differential equations for the probability distribution. The system 
of differential equations take the form of master equations for a jump Markov process and 
Fokker-Planck equations (FPE) for a continuous Markov process. For the way in which 
this happens, the reader is referred to [? ] and [? ? ]. For a Markov process that is 
made up of both jump and continuous parts, the differential equation takes the form of 
the differential Chapman-Kolmogorov equation (dCKE) which has been derived in [? ]. 
The derivation is involved and requires the introduction of an arbitrary function, which 
leads to boundary restrictions on the probability distribution. As part of this review, we 
present a novel and more concise derivation of the dCKE. Since most of the mathematical 
foundations for stochastic models have been developed by physicists and mathematicians, 
we hope that our derivation makes the theory more accessible to the uninitiated researcher 
in the field of systems biology. We choose Markov processes as a framework, since more 
realistic approaches for modelling intracellular processes must take into account factors 
such as heterogeneity of the environment, macromolecular crowding [? ? ] and anomalous 
diffusion [???], to name a few. Anomalous diffusion is described by fractional Fokker- 
Planck equations [? ]. Such treatments require advanced mathematical formalisms which 
are beyond the level assumed in this paper. 

The focus of the present paper is neither a comprehensive review of stochastic approaches 
(See [? ] for a recent survey, [? ? ] for a recent theoretical analysis, [? ] for a discussion 
of the role of stochasticity in cell biology) nor a comparison of the two approaches (e.g. [? 
? ? ]). Instead, we review the formal relationships between the equations referred to in 
the systems biology literature. We thereby, try to establish a "family tree" for stochastic 
models in systems biology, providing explicit derivations of their formal relationship and 
clarifying assumptions involved in a common framework (See Figure [2]). In the following 
section we focus on the origin of the chemical master equation CME (a special form of 
the master equation for systems governed by chemical reactions) within the framework 
of Markov processes. Such generalisation provides a clearer picture of how the various 
stochastic approaches used in systems biology are related within a common framework. 

2 Markov processes 

Markov processes form the basis for the vast majority of stochastic models of dynamical 
systems. The three books by ? ], ? ] and ? ] have become standard references for the 
application of Markov processes to biological and biochemical systems. At the centre of 
a stochastic analysis is the so-called Chapman-Kolmogorov equation (CKE) that describes 
the evolution of a Markov process over time. From the CKE stem three equations of 
practical importance: the master equation for jump-Markov processes, the Fokker-Planck 
equation (FPE) for continuous Markov processes and the differential Chapman-Kolmogorov 
equation (dCKE) for processes made up of both the continuous and jump parts. A nice 
mathematical (but non-biological) account of these equations can also be found in [? ]. 
Gardiner, van Kampen and Gillespie take different approaches to derive these equations: 

• Gillespie derives the FPE and the master equation independently from the CKE and 
for the one-dimensional case only in [? ]. In [? ? ] he extends the derivations to 
multidimensional cases. 

• ? ] derives the master equation from the CKE for a one-dimensional Markov process. 
The FPE is given as an approximation of the master equation by approximating a 
jump processes with a continuous one. The same approach is adopted in [? ]. 
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However, this should not mislead the reader to conclude that the FPE arises this 
way. In fact, FPE is defined for a continuous Markov process. 



• ? ] derives first the dCKE from the CKE for a multidimensional Markov process 
whose probability distribution is assumed to be contained in a closed surface. The 
FPE and the master equation are given as special cases of the dCKE. 

We start with a review of basic concepts from probability theory required to read our proof. 
This is followed by a brief derivation of the CKE and its graphical interpretation. From 
the CKE we derive the dCKE and interpret its terms to show how the FPE and the master 
equation appear as special cases of the dCKE. Finally we show that the CME is just a 
special form of the master equation for jump processes governed by chemical reactions. 

A random variable X describes a random event by assigning it values x (called states) 
from a set S (called state-space) and defines a probability distribution over this set. The 
set S may be discrete, continuous or both. The probability distribution, in the case of a 
discrete state-space S = {n}, is given by a set of probabilities p n such that 

Prob {X = n} = p n . 

In the case of a continuous state-space S = {x}, the probability distribution is given by a 
non-negative function p(x) such that 

Prob {x < X < x + dx} = p{x)dx . 

In the literature, p n is referred to as a probability mass function (p.m.f) and p(x) as a 
probability density function (p.d.f). The delta function 5(x) defined by 

J ' dxf(x)5{x-c) = /(c) (1) 
s 

for any function / of x, allows us to write the discrete distribution as a special case of the 
continuous case. Specifically we write 

P( x ) = ^2PnS(x - n) 



and note that 
Prob 



{X = n} = / dxp(x) = / dx S ^p m 5(x — m) = p n 



x<n<x+dx x<n<x+dx 



which is what we would expect in the discrete distribution. Since the discrete distribution 
can always be derived easily from a continuous one, we use hereafter the latter. When 
dealing with dynamical systems, the probability distribution evolves over time. This leads 
to the notion of a stochastic process, that is, a system in which random variables are 
functions of time, written as X(t). The states are scalars for a one-dimensional system 
and vectors for a multidimensional system. Note that, while the derivations of the master 
equation and the Fokker-Planck equation in [? ? ] are for the one-dimensional case only, 
we here present a general treatment and derive all our results for multidimensional systems. 
The probability distribution for an iV-dimensional stochastic process 

Mt) = (x 1 (t),...,x N (t)) 
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N >| 

Pi Xi < Xi(t) < Xi + dxi > 
i=l J 



. . . ,X]y,t)dxi ■ ■ ■ dxN ■ 



, 1= 

To simplify the notation, we use a short form 



Prob {x < X(t) < x + dx} = p(x, t)dx . 

More useful will be the conditional probability density, p(x, t |x',t'), defined such that 

Prob {x < X(t) < x + dx | X(t') = x } = p(x, t | x , t')dx . 

When t > t', p(x, t 1 x', t')<ix is called transition probability. Since it is much easier to work 
with densities p(-) rather than probabilities Prob{-}, we shall use densities p(-), but abuse 
the terminology by referring to it as "probabilities". 

Essentially a Markov process is a stochastic process with a short term memory. Math- 
ematically it means that the conditional probability of a state is determined entirely 
by the knowledge of the most recent state. Specifically for any three successive times 
to < t < t + At , one has 

p(x, t + At | z, t; xo, to) = p(x, t + At | z, t) 

where the conditional probability of x at t + At is uniquely determined by the most recent 
state z at t and is not affected by any knowledge of the initial state xo at to- This Markov 
property is assumed to hold true for any number of successive time intervals. To see how 
powerful this property is, let us consider the factorisation of the joint probability 

p(x, t + At ; z, t) = p(x, t + At | z, t)p(z, t) . 

Making both sides conditional on (xo,to) will modify this equation as 

p(x, t + At ; z, 1 1 x , t ) = p(x, t + At \ z, t; x , t )p(z, t \ x , t ) 

which, by using the Markov property, reduces to 

p(x, t + At ; z, 1 1 x , to) = t + At \ z, t)p(z, t \ x , t ) . (2) 

The last equation shows that the joint probability can be expressed in terms of transition 
probabilities. Recall the following rule for joint probabilities 

p(x) = J dyp(x,y) (3) 

which says that summing a joint probability over all values of one of the variables eliminates 
that variable. Now integrating j2|) over z and using ([3|), we arrive at the so-called Chapman- 
Kolmogorov equation (CKE) [? ]: 



p(x, t + At | Xq, to) = J dzp(x, t + At|z, t)p(z, t|xo, to) • 



(4) 



This equation expresses the probability of a transition (xo — > x) as the summation of 
probabilities of all transitions (xo — ► z — > x) via intermediate states z. Figure Q] illustrates 
the basic notion of a Markov process for which the CKE provides the stochastic formalism. 
When the initial condition (xo, to) is fixed, which is assumed here, the transition probability 
conditioned on (xo,to) is the same as the state probability: 

p(x,t) =p(x,t|x ,t ). 
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3 Derivation of the dCKE 



The CKE serves as a description of a general Markov process, but cannot be used to deter- 
mine the temporal evolution of the probability. Here we derive from the CKE a differential 
equation which will be more useful in terms of describing the dynamics of the stochastic 
process. Referred to as the differential Chapman-Kolmogorov equation (dCKE) by ? ], 
this equation contains the CME as a special case. This derivation is for a multidimensional 
Markov process. We start with the definition of a time-derivative for a probability density 
but place no restrictions on the probability distribution. ? ] instead starts with the expec- 
tation of an arbitrary function which results in integration by parts and consequently the 
need to assume that the probability density vanishes on the surface of a region to which 
the process is confined. We do not need such as assumption because of the simplicity of 
our approach. The master equation gives the jump part of the Markov process while the 
Fokker-Planck equation gives the continuous part. 

Consider the time-derivative of the transition probability 

-p(x, t) = ttm — {p(x, t + At) - p(x, t) } , (5) 

where differentiability of the transition probability with respect to time is assumed. Em- 
ploying the CKE (J4j) and the normalisation condition 



J dzp(z,t + Ai|x,t) = 1, 

s 



since p(z, t + At | x, t) is a probability, ([5]) can be rewritten as 



d_ 
dt 



p(x, t) = j dz |p(x, t + At | z, t)p(z, t) — p(z, t + At | x, t)p(x, t) | 



Let us divide the region S of integration into two regions based on an arbitrarily small 
parameter e > 0. The first region ||x — z|| < e corresponds to a continuous state process 
and the above derivative in this region will be denote by I\. Here ||-|| denotes a suitable 
vector-norm. The second region ||x — z|| > e corresponds to a jump process and the above 
derivative in this region will be denote by I2. Since (|6|) gives the derivative in the whole 
region S, we can write 

|U(x,t) = J 1 + I 2 , (6) 



dt' 



where 



h = ^i m Q ^ j dz jpCM + At \z,t)p(z,t) -p{z,t + At I x, t)p(x, t) j . 



||x— z||<e 

and 



I2 = Jim J dz I p(x, t + At \ z, t)p(z, t) — p(z, t + At | x, t)p(x, t) > . 

||x-z|]>e ^ ' 

In the first region ||x — z|| < e, the integrand of I\ can be expanded in powers of x — z 
using a Taylor expansion. Setting x — z = r, we can write, 

Ii=2im o ^ j dr |p(x,t + At | x - r, t)p(x - r, t) -p(x-r,t + At |x,t)p(x,t)j . (7) 
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In order to expand the integrand more easily into a Taylor series, let us define a function 

/(x; r) 4 p( x + r, t + At | x, t)p(x, i) 

so that the integrand in ([7j) becomes /(x— r; r) — /(x; — r), which, after a Taylor expansion, 
becomes 

(9/(x'r) 1 3 2 /(x - r) 
-/(x; -r) + /(x;r) + }(-n) — ^ h „ T^f"^- — — ^ ^ higher-order terms . 

1 i, j J 

The integrals of the first two terms cancel because of the symmetry J /(x; r)dr = J* /(x; — r)dr, 
when the integral is over all the positive and negative values of r in the region. Thus, we 
have 



/1 = lim , 

At^O At 



_d_ 

1 dxi 



p(x + r, t + At I x, t)p(x, t) 



\\r\\<e 

1 

2 



+ 9 Yj nT 3 



2 



dxidxj 



p(x + r, t + At I x, t)p(x, t) + higher-order terms 



For the state increments X{(t + At) — -Xj(t), recognising the (conditional) expectations 
(X;(t + At) - Xi(t) I X(t) = x) = J drrip(x + r,t + At|x,t) 



and 



([X;(t + At) - Xi{t)} [Xj(t + At) - Xj(t)} I X(t) = x) 



J dr TiTj p(x + r, t + At | x, t), 



r <e 



we refer to the differentiability conditions for continuous processes, i.e., ||x — z|| < e [? , 
section 3.4]: 

{X i (t + At)-X t (t)\X(t)=x) 



lim 

At-^0 



At 



^(x,t) + o(e) (8) 



Hm {[ X l{ t + At)-X M[ X j{ t + At)-X M |X(t)=x) = ^ (x>f)+o(e) (9) 

_A / *0 _A / 

where o(e) represents vanishing terms, such that lim e ^oo(e)/e = 0. The higher-order 
terms involve higher-order coefficients which must vanish. To see that, for the third-order 
coefficient 



lim 

At 



^^t / dvriTjrkpi-x. + r, t + At |x,t) = C ijk (x,t) + o(e) . 



k <e 



However 



lim 

At^O 



J drrirjr k p(-x + r,t + At|x,t) 
1 f 

At / ^ r r * r i^( x + r ' * + I x > ^) 



k <e 



< ||r|| lim 
At 

lkll< e 

< e[fl -(x,t)+o(e)] 

< o(c). 
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Hence C(x, t) must vanish. The vanishing of higher-order coefficients follows immediately. 
In physics, the vector A(x,t) is known as the 11 drift vector" and the matrix i?(x, t) as the 
" diffusion matrix". This terminology is suggested by the observation that, given X{t) = x, 
the state increment vector X(t+dt) — X(t) for a continuous process has a mean approaching 
A(x,t)dt and a covariance approaching B(x,t)dt, as e approaches zero. This also suggests 
the following update rule, for e — > and under assumptions given in [? , section 3.5.2], 



X(t + dt) = X(t) + A(X(t),t)dt + [B(X(t),t)dt] 



1/2 



(10) 



which is a form of the "Langevin equation" [? ]. We remark here that ([8]) and (jQj) are 
postulated here for mathematical convenience. A more rigorous justification is given in [? 
]. Subject to the differentiability conditions ((HI) and ((9|), we see that as e — > 0, 



o 2 



dxi . 

Next we work out the jump probability rate 
h 



dxjdxn 



By-(X 3 t)p(x,t) 



[ID 



hJ 



lim j dz ^p(x, t + At | z, t)p(z, t) — p(z, t + At | x, t)p(x, t) ^ 



||x— z || >e 

We will use the differentiability condition for jump processes, i.e., ||x — z|| > e [? , section 
3.4]: 

lim -r-p(x, t + At I z, t) = W(x I z, t) , 
At^o At v 1 ' 

where W(x|z,t) is called the transition rate for the jump (z — > x). Subject to this 
condition, we see that as e — > 0, the region of integration approaches 5, leading to 



/ 2 



/ 



dz 



W (x | z, t)p(z, t) - VF(z | x, t)p(x, t) 



[12) 



Adding (fTH and (TT2l) . we can rewrite J6]) to arrive at the dCKE: 



d_ 

dt 



p(x, t) 



a 2 



2 ^ dxidxj - 



Sy(x,t)p(X,t) 



+ y dz [w(x | z, t)p(z, t) - W(z | x, t)p(x, t) 



(13) 



We now have a differential equation characterising the dynamics of the probability distribu- 
tion p(x, t), that is the probability of a state at any time, starting from a given initial proba- 
bility distribution. This completes our derivation of the differential Chapman-Kolmogorov 
equation. Differences between our derivation and those available in the literature are de- 
scribed in Section [7] (Conclusions). The following section will classify Markov processes 
based on this dCKE. This is followed by a derivation of the chemical master equation. 
Finally, in Section [61 we discuss the use of the master equation in systems biology. 



4 Classification of Markov processes based on the dCKE 

Being a linear differential equation, the dCKE is more convenient for mathematical treat- 
ment than the original CKE. More importantly, it has a more direct physical interpretation. 
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The coefficients A(x,t), B(x,t) and W(z|x, t) are specified by the system under consid- 
eration, and thus the solution of the dCKE gives the probability distribution for the state 
of the given system [? ]. The original CKE, on the other hand, has no specific informa- 
tion about any particular Markov process. We now interpret the different terms of (fT3l) , 
Following [? ? ] we first consider the case 



reducing the dCKE to 



Bij(x,t) = W(x|z,t) = W(z|x,t) = 0, 



which is a special case of the so-called "Liouville equation", describing a deterministic 
motion (See [? , section 3.5.3]): 

j*(t) = A(x,t). 

This is the simplest example of a Markov process. 
Next, if A(x,t) = B(x,t) = 0, the CKE reduces to 



d f r 

— p(x, t) = / dz [W{x | z, t)p(z, t) - W{z | x, t)p{x, t) 



(14) 



This is called the "master equation" describing jump-Markov process with discontinuous 
sample paths. 

Next, if W(x | z, t) = W(z | x, t) = 0, the CKE reduces to 



J^(x, t) = - ^7 h( x > *)^( x ' *)1 + \ E 



c) 2 



dxidxj 



Bij(x,t)p(x,t) 



which is called the "Fokker- Planck equation" (FPE) and is equivalent to the Langevin's 
equation (fTUI) under the conditions given in [? ? ? ]. The corresponding process is 
known as a diffusion process which is continuous but not deterministic. This shows that 
the FPE is originally defined for a continuous process. However the FPE can also arise 
as an approximation of the master equation when the jumps of the corresponding discrete 
process are assumed to be small [? ? ]. 

Finally we consider the case where the diffusion matrix B(x,t) = 0, which leads us to 



d d r 

— p(x, i) = - £ — Ai(x, t)p(x, t) 



dxi 



+ 



dz 



VF(x | z, t)p(z, i) - W(z | x, i)p(x, i) 



which is called the "Liouville master equation" (LME) in [? , chap. 1] and describes 
a piecewise deterministic process with sample paths consisting of smooth deterministic 
pieces interrupted by instantaneous jumps. One way in which the LME arises is when an 
originally jump Markov process is approximated by a process with discrete and continuous 
components [? ? ]. 

In the most general case, when none of the quantities A(x.,t), B(x,t) and W(z\x,t) 
vanish, the dCKE may describe a process whose sample paths are piecewise continuous, 
made up of pieces which correspond to a diffusion process with a nonzero drift, onto which 
is superimposed a fluctuating part. 
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5 The chemical master equation 



Consider a Markov process with a discrete state-space. The master equation for this 
discrete process can be obtained from ([HJ to give 

m 

where m the intermediate state, and n the final state. Since p n (t) is a probability (and 
not a density), the integral j s is replaced with the summation X^m- We can rewrite this 
equation in terms of jumps r = n — m, 

J^n(t) = Y [ W ( n I n " r > *)Pn-r(*) W(n + r | n, i)p n (i)] , (15) 

r 

where we have used the symmetry ^2 r 4>{— r) = X^r^( r )' f° r an arbitrary function 
when writing the second summand. Now consider an iV-component and M-reaction bio- 
chemical system. Let i label the different components (chemical species) and j label dif- 
ferent reaction channels. The copy number of ith component at the variable time t will be 
denoted by X,i{t) which takes values rii from the set of whole numbers. Each occurrence of 
j-th reaction channel changes the copy number rij of ith component by an amount Uij. The 
elements z/jj form the stochiometric matrix v whose jth column will be denoted by Uj. It is 
assumed that the species are distributed homogeneously (well mixed) in a closed system of 
constant volume O at a constant temperature. This essentially assumes that changes only 
depend on the current state (Markov property) and that we can avoid spatial considerations 
[? ? ? ] and macromolecular crowding [? ]. However, since diffusion may not always be 
rapid, spatial considerations become important when dealing with intracellular processes 
[???]. Here we are interested in a stochastic formulation which dates back to the 
initial work by ? ]. Under the stated assumptions, the vector X(i) = (Xi(t), . . . ,Xjf(t)) 
taking values n = (m, . . . , tin) is a continuous time Markov process. The jump sizes v are 
determined by the stoichiometry and molecularity of the reactions and, therefore, can only 
take values from the set {ui, . . . , um} of the elementary changes. Thus, for our system of 
chemical reactions, (fT5l) becomes 



d M r 

3=1 

Since Uj is uniquely defined for a reaction Rj, we introduce a simpler notation 

Oj(n) = W(n + uj | n, t) 
to rewrite the above master equation as: 

M 

dt 



d 

JT,Pn{-t) = Yl M n ~ V j)Pn-Vi (*) ~ a j( n )Pn(t)] , (16) 



3=1 

which is referred as the chemical master equation in the systems biology literature [???]. 
This shows that the CME is just a special form of the master equation for jump processes 
governed by chemical reactions. The coefficient a,j is referred to as the reaction propensity 
and is interpreted such that a,j(n)dt gives the probability of jth reaction occurring in the 
time interval [t, t+dt) from state n at time t. In the stochastic setting of ? ], the jth reaction 
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channel is characterised by a stochastic rate constant Cj such that Cjdt gives the probability 
that a particular combination of molecules will react according the j'th channel in the next 
infinitesimal interval [t,t + dt]. The propensity aj(n) is thus Cj times the number hj(n) of 
different possible ways in which molecules can combine to react according the j'th channel. 
Since this equation is difficult to solve analytically or even numerically, several attempts 
have been made to avoid a direct solution or simulation of the CME. The most successful 
implementation is the stochastic simulation algorithm (SSA) which originated from work 
by ? ] but it was Gillespie who pioneered its use for generating sample paths of chemical 
reaction networks [? ? ]. The SSA is widely used in systems biology [?????]. Figure 
[2] provides an overview of stochastic models and interrelationships referred to here. 

This completes our formal analysis and we now return to application of the CME in 
systems biology. 

6 Master equations in systems biology 

The chemical master equation (fT6j) is the basis for most stochastic models in systems biol- 
ogy. For complex systems, involving large numbers of chemical species and reactions, the 
computational cost may be considerable. For this reason modifications to the algorithm and 
strategies to simplify the model prior to a computer simulation have been suggested. The 
efforts to reduce the computational complexity of stochastic approaches can be grouped 
into approximate stochastic methods and hybrid methods. Approximate stochastic meth- 
ods try to speed up the simulation by compromising exactness of the master equation 
whereas the hybrid methods treat parts of the system deterministically and other parts 
stochastically. What follows is a brief discussion of the systems biology literature, using the 
CME and SSA, including strategies that have been developed to reduce the computational 
complexity. 

6.1 Approximate stochastic methods 

In [? ] Gillespie presents an approximate and thereby faster simulation method known as 
the r-leap method. Instead of simulating individual reactions, the number of reactions of 
each type in a sequence of short time intervals is simulated. Optimal ways to select the 
leap-length of the intervals have been investigated in [? ? ]. In [? ] Rao et al. propose a re- 
duction of a chemical system by partitioning molecular species into slow (primary) and fast 
(intermediate) molecular species. Assuming that the intermediate species (conditional on 
the primary species) are Markovian, they apply the quasi-steady-state assumption (QSSA) 
which essentially assumes that the conditional probability distribution of the intermedi- 
ate species is time-invariant (i.e., it has reached a steady-state); thereby eliminating these 
species from the chemical master equation. A modified version of Gillespie algorithm is 
subsequently proposed to simulate the resulting reduced system for slow species. In [? 
? ] Haseltine and coworker also use the partitioning method for reduction, but partition 
chemical reactions into fast and slow reactions. The fast reactions are approximated by 
using Langevin or deterministic equations. The coefficients of the reduced chemical master 
equation are influenced by the fast reactions. The authors propose simulation algorithms 
for the slow reactions, subject to constraints imposed by fast reactions. The idea of parti- 
tioning to speed up slow-scale simulation has also been used in [????].? ] present a 
probabilistic steady-state approximation that separates the time scales of an arbitrary re- 
action network, detects the convergence of a marginal distribution to a quasi-steady-state, 
directly samples the underlying distribution, and uses those samples to predict the state 
of the system. ? ] propose that, in case of higher dimensions, the master equation could 
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be approximated by the FPE and then discretised in space and time by a finite difference 
method. They demonstrate the method for a four- dimensional problem in the regulation 
of cell processes and compare it to the Monte Carlo method of Gillespie. ? ] use the CME 
to analyse a negative feedback system composed of two species regulating the synthesis 
of each other. In [? ] Paulsson uses a variant of the fluctuation-dissipation theorem to 
give a generic expression for noise arising from different cellular processes, applying the 
theorem to a simple generic model representing simple gene expression. Paulsson uses the 
notion of an f2-expansion of the master equation, i.e., a Taylor series expansion in powers 
of where Q is a system size parameter [? ]. The first- and second- order terms of the 
expansion reproduce the macroscopic rate equations and realise the fluctuation-dissipation 
theorem respectively. An equation that decomposes the intrinsic and extrinsic noise con- 
tributions is derived to simplify the analysis. ? ] give a general method to simplify the 
master equation in a linear noise approximation (LNA), obtained through an O-expansion 
of the master equation. They derive the LNA for the stationary state of a general system 
of chemical reactions and use it to estimate sizes, correlations and time scales of stochastic 
fluctuations. They demonstrate that the LNA allows a rapid characterisation of stochastic 
properties for intracellular networks over a large parameter space. They also show that 
the LNA can be made more accurate in cases where fast variables can be eliminated from 
the system. In [? ] the same authors use the LNA of the master equation to characterise 
intracellular metabolite fluctuations. In both publications, the results of the LNA are 
compared to simulation of the master equation. 

6.2 Hybrid methods 

? ] present a simulation approach for hybrid stochastic and deterministic reaction models. 
The system is adaptively partitioned into deterministic and stochastic parts based on a 
given criteria at each time step of the simulation. They present two algorithmic schemes. 
The 'direct hybrid method' explicitly calculates which reaction occurs and when it occurs. 
The 'first and next reaction method' generates a putative time for each reaction. The reac- 
tion corresponding to the smallest time is chosen to occur; the state according updated and 
the process repeated. ? ] derive the first-order partial differential equations for probability 
distribution function from stochastic differential equations describing approximate kinetics 
of a single cell. Resulting equations are used to calculate mRNA-protein distribution in the 
case of single gene regulation and the protein-protein distribution in the case of two-gene 
regulatory systems. In [? ] the same authors present a hybrid stochastic and determin- 
istic treatment of the NF-ftB regulatory module to analyse a single cell regulation. They 
combine ordinary differential equations, used for the description of fast reaction channels 
of processes involving a large number of molecules, with a stochastic switch to account for 
the activity of the genes involved. ? ] makes two approximations to the exact stochastic 
description of a gene regulatory network: the continuous approximation considering only 
the stochasticity due to the gene activity; and the mixed approximation attributing the 
additional stochasticity to the mRNA transcription/decay process. The underlying dis- 
tribution is then described by a system of partial differential equations derived from the 
dCKE after specific assumptions on the coefficients A, B and W. ? ] take a different 
approach by deciding on the fly which approach to use. Studying the stochastic simulation 
of signal transduction via calcium, the authors observe that the transition from stochastic 
to deterministic occurs within a range of particle numbers that depends the phase space 
of the system. 

Yet another approach is given in [? ] where a novel multidimensional stochastic frame- 
work is proposed to model multi-gene expression dynamics. Inspired by the dCKE, they 
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propose a new experimental scheme which will measure the instantaneous transition prob- 
abilities. Given experimental data, one could obtain the coefficients of the dCKE, which 
can then be solved to obtain the distribution and the corresponding moments and corre- 
lations. ? ] analyse noise in a negatively feedback-regulated transcription factor (TF) 
and the effects of the feedback loop on a gene repressed by the same TF using a modified 
Gillespie algorithm (provided by SmartCell). The authors find that within a certain range 
of repression strength, the negative feedback loop minimises the noise whereas outside this 
range, noise is increased. It is proposed that this may arise from plasmid fluctuations. 

7 Conclusions 

The application of stochastic models is usually motivated by uncertainty arising from vari- 
ability. In the engineering and physical sciences the variability arises mainly from measure- 
ments. In systems biology the variability arises mainly from the complexity of intracellular 
processes. By complexity we mean the fact that in a particular biological network we are 
forced to eliminate many variables that are influencing the observation we make. While 
the assumptions of constant temperature, pH level, volume and water balance may not 
worry most modellers, the large number of unmodelled variables may be of greater con- 
cern. Implicit in most stochastic models is the assumption of a well-mixed homogeneous 
environment in which there are more non-reactive collisions than the reactive ones. More 
realistic spatial representations are therefore an increasingly important research theme in 
systems biology. In gene expression a very small number of molecules controls potentially 
very large molecular populations, suggesting hybrid approaches to combine stochastic and 
deterministic formalisms. Many cellular processes, for instance the differentiation of stem 
cells, are multistable systems for which state-of-the-art single-cell measurements are pro- 
viding increasingly valuable data with nanometre and millisecond resolution. The advance 
of these technologies allows us to monitor the transcription of individual genes, leading 
also to a demand for advanced stochastic modelling and simulation. 

Motivated by applications of stochastic models in systems biology, we described a prob- 
abilistic framework based on Markov processes to represent biochemical reaction networks. 
We provided a novel derivation of the differential Chapman-Kolmogorov equation for a 
general Markov process made up of both continuous and jump processes. Then we re- 
viewed the formal relationships between the equations referred to in the systems biology 
literature, establishing a "family tree" for stochastic models in systems biology, providing 
explicit derivations of their formal relationship and clarifying assumptions involved in a 
common framework (See Figure [2]). Our derivation starts with the definition of a time- 
derivative unlike Gardiner who starts with the expectation of an arbitrary function. We 
place no restrictions on the probability distribution, whereas Gardiner assumes it to be 
confined to a region that has a surface, and the probability being zero on the surface. 
The master equation gives the jump part of the Markov process while the FPE gives the 
continuous part. The derivation of FPE in ? ] and ? ] involves approximation of the 
master equation by investigating the limit of a jump process. 
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Figure 1: Graphical interpretation of the Chapman-Kolmogorov equation. The probabil- 
ity of a transition (xo — > x) can be obtained by summing probabilities of all 
transitions (xq — > z — > x), via intermediate states z. 
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Figure 2: Interrelationships for the stochastic models and their simulation which are cov- 
ered in this paper. The coefficients A, B, W respectively refer to the drift-vector, 
diffusion-matrix and the transition-rate in the dCKE. QSSA stands for Quasi- 
Steady-State Assumption. 
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